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Abstract. While there has been a large interest in studying the role of 
dissolution-driven free convection in the context of geological sequestration, 
the contribution of forced convection has been largely ignored. This manu- 
script considers CO2 sequestration in saline aquifers with natural background 
flow and uses theoretical arguments to compute the critical background veloc- 
ity needed to establish the forced convective regime. The theoretical arguments 
are supported by two dimensional high-resolution numerical simulations which 
demonstrate the importance of forced convection in enhancing dissolution in 
aquifers characterised by low Rayleigh numbers. 



1. Introduction 

Geological storage of carbon dioxide in deep saline aquifers is considered to be one 
of the key technological solutions to mitigating carbon emissions in the near future 



(IPCC 2005 ). Under the thermodynamic conditions at which CO2 is injected into 
subsurface formations, CO2 exists in a supercritical state with a density that is 
lower than that of brine in the medium. Due to this difference in density, CO2 
rises to the top of the domain and containment is achieved by the presence of an 
overlying low-permeability caprock. 

Due to the tremendous importance of understanding the processes leading to 
storage security of injected carbon dioxide, there has been a large body of published 
literature exploring the different trapping mechanisms. The physical processes that 
play a role at increasing timescales are (i) physical trapping due to the caprock, (ii) 
residual trapping of bubbles of CO2 in the porous rock, (Hi) dissolution trapping 
due to the slow dissolution of CO2 into the brine and (iv) mineral trapping due 
to chemical reactions between dissoved carbon dioxide and host rocks. The work 
presented in this manuscript is concerned with dissolution trapping mechanism 
under different aquifer flow conditions. Experimental data has shown that the 
dissolution of CO2 into brine produces a weak density contrast, with CO2 saturated 
brine being slightly denser than unsaturated brine. Most recent work in this area 
has focused on the free convective mixing that can arise due to this difference in 
density. 

One of the effects that has so far been neglected is the contribution of background 
flow to the rate of dissolution of CO2 into the brine in the aquifer. The only 



exception we are aware of is the work by Hassanzadeh et al. 2009 , where the 
authors have analysed the role of additional dispersion due to the background flow. 
For a finite sized CO2 plume in the supercritical state, the background flow leads 
to the development of a boundary layer in which horizontal advection balances the 
diffusion of dissolved carbon dioxide. We will henceforth refer to the dissolution 
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occuring due to the background flow as forced convection. A compilation of reservoir 



characteristics for potential sites in western Canada by Bachu et al. 2004 suggests 
that the Rayleigh numbers are going to be fairly low for a large number of storage 
sites. Recent computations (Rapaka et al. 2008j) as well as experiments (Backhaus 



et al. 2011 ) in analogous systems have shown that the free convective mechanism 
(fingering) can take very long to establish under such conditions. It is possible 
that forced convection can play the dominant role in enhancing storage security 
in this scanario. However, in the context of geological storage of carbon dioxide, 
the conditions under which forced convection is the preferred state have not yet 
been explored. In this manuscript, we will provide simple theoretical arguments 
for the background flow conditions under which forced convection establishes itself 
and support these predictions with high resolution numerical simulations. 

2. Governing Equations and Solution Procedure 

We consider the dynamics of CO2 dissolution and transport in a aquifer of depth 
H* with an overlying source of CO2 of length L* . We assume the presence of a 
horizontal background aquifer flow Uq = (ug , 0) , the magnitude of which determines 
whether the system is in a free or forced convective regime. 

The equations needed to describe the system are those describing conservation of 
mass, momentum (Darcy's law) and an advection-diffusion equation for the concen- 
tration of dissolved CO2. Experimental data has shown that the density of brine 
increases linearly with the concentration of dissolved CO2, with a maximum in- 
crease (under saturated conditions) of the order of 1% (Pruess and Spycher 2006| ). 



d{(t>C* 
dt* 



Hence, we utilize the Boussinesq model, where the density of brine is taken to be a 
constant everywhere except in the buoyancy term. The governing equations are: 

!^u*^~V*P*+p}gey, (1) 
V* • u* = 0, (2) 

V* ■ {u*C*) =V* ■ (t)D*V*C*, (3) 
p}=p*o + Ap* (^^) , (4) 

where, u* = {u* , v*) is the velocity of the fluid, C* is the concentration of dissolved 
carbon dioxide, P* is the fluid pressure, g is the acceleration due to gravity, k* is the 
permeability of the rock, p^ is the density of unsaturated brine, Ap* is the density 
difference between brine saturated with CO2 and pure brine, p* is the viscosity of 
brine, (jj is the porosity of the medium and D* is the diffusivity of CO2 in brine. 

We now non-dimensionalize the equations choosing a convective velocity scale 
Ur — k* Ap* g / p.* , the system depth as a length scale Lr — H* , and a convective 
time scale = Lj./Ur- Concentrations are non-dimensionalized using Cq, the 
maximum solubility of CO2 in brine at the conditions in the aquifer. The non- 
dimensional equations obtained are: 

u = -VP + Cey (5) 
V • u = (6) 

^ + v-(uc).±v>a (7) 
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Figure 1 . Schematic illustration of the problem geometry consid- 
ered in this manuscript. The red part of the top boundary denotes 
the footprint of CO2 plume in supercritical state. A mean hor- 
izontal flow of magnitude uq is taken to be flowing through the 
system. 



Ra = 



(8) 



The non-dimensional paramerer Ra is the ratio of the convective flux to the diffusive 
one: 

k*Ap*gH* 

The initial conditions are given by C ~ Q,u — uq and v — Q everywhere. The 
CO2 concentration is fixed on the top boundary to be C = 1 for a source region 
of length L, and a no-flux boundary condition is applied for the remainder of the 
top boundary and the entire bottom boudary. On the left boundary, we have 
C ^ 0,u ~ uq and v — and an outflow boundary condition is applied on the right 
boundary. A schematic of the geometry considered is shown in Fig. ([I]). 

Depending on the magnitude of the background flow uq, CO2 dissolution can 
be enhanced by free convection (low uq, sufficiently high Ra), forced convection 
(high Mo) or be a purely diffusion-driven process (low uq, low Ra). The case of free 
convection (especially in the absence of any background ffow) has received a great 
amount of attention in the recent past (see for instance Ennis-King et al. [2005 , Xu 

, [Hassanzadeh et al. 2007] , Rapaka et al. [2008] and 
) . There is an extensive body of published research on 
forced convection in porous media over the last four decades, and it is not possible 
to present present an overview of even the most important contributions within 



et al. 


2006 , 


Riaz et al. 


2006 


Slim and Ramakrishnan 


2010 



the scope of this paper. The reader is referred to the book by |Nield and Bejan 



2006 



where a comprehensive account of the existing body of work is presented. 
However, in the context of geological storage of carbon dioxide, we are not aware 
of any manuscript that points out the conditions under which forced convection 
plays a dominant role in enhancing the dissolution of injected CO2. For sufficiently 
large background flow rates, we expect to see a forced convective regime achieving 
a steady-state dissolution rate. Using the commonly made assumption of << 
in the boundary layer and assuming that the basic velocity fleld is of the form 
u — uq,v — 0, the steady state advection diffusion equation can be reduced to: 

dC 1 d^C 
'''d^ = R-aW 

Looking for similarity solutions of the form 1 — C{x,y) = /(^) in terms of the 
similarity variable ^ ~ y/^/x, it is easy to show that /(^) satisfles the following 
differential equation: 

1 



/" + -Rauof^ = 



(10) 
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where, primes denote differentiation witfi respect to the similarity variable. The 
solution of this equation can be shown to be: 



fiO = erf{^^Rauo/2) 
giving the dissolution rate at the top boundary as 



dC 
dy 



iy = o) 



^/RauQ 



(11) 



(12) 



to L, we get the total dissolution rate to be 
= UqL* /(pD is the Peclet number based on the 



Integrating this expression from x = 
Nu = -i/RauQL — %/ Pe, where Pe 
size of the CO2 plume. 

It can be seen that the total dissolution rate under conditions of forced convection 
scales with the length of the plume as VL. In the free convection regime however, 
the total dissolution rate scales linearly with L. We expect the presence of a critical 
velocity Uc{L) such that, for slow background flows characterized by Uq < Uc{L), 
the system will be in a state of free (mixed) convection, whereas for uq > Uc{L), 
we will transition to the forced convection regime. 

To determine the critical background flow velocity Uc{L), we need an expression 
for the mean dissolution rate as a function of Rayleigh number under free convective 
conditions. Recent studies have looked at the time-averaged Nusselt number under 
free convective conditions and proposed scalings of the form Nu ~ Ra". |Neufeld| 
have used scaling arguments and experiments to suggest a = 0.8, 
have used experiments to obtain a best-fit exponent 



et al. 



2010 



Backhaus et al. 2011 



while 

of a = 0.76. However, the experiments by both authors were performed under 
conditions of very high Rayleigh numbers, of the order of lO'* or higher. In this 
study, we will use a simpler model given by 

RaL 



40 



(13) 



which is the well-known relation between Nusselt number and Rayleigh numbers for 
steady convection when Ra is close to the critical Rayleigh number (Elder 1967| ). 
This relation is known to overpredict the Nusselt number for large Ra, but is in 
good agreement with experimental data for Ra < 250. Since the focus of this 
manuscript is on aquifers with lower permeabilities, we use this relation as a useful 
approximation whose accuracy will be examined a posteriori. Requiring the free 
and forced convective fluxes to agree at uq — Uc{L) , we get 



RaL 
1600' 



(14) 



We now look at some values for the parameters characteristic of real formations. 
We use data from Bachu et al. 2004 , where a vast amount of information about 



the potential s ites in western Canada is presented. Using the data from |Pruess and| 
Spycher 2006 , we estimate the density increase Ap* to be 5kg /m^ (corresponding 
to a salt mass fraction of 0.1). Assuming k* = IQmD, g — lOm/s^, H* — 5Qm, /i* — 
0.5mPa.s, (f) = 0.2 and D* = 1 x 10^^, we get the Rayleigh number to be Ra = 250. 
The convective velocity scale is Ur — k*Ap*g/fi* — icm/year. Assuming a CO2 
plume size of 2km, we can compute the critical background flow velocity needed 
to transition to a forced convective regime as u* = 20cm /year (where the asterisk 
denotes dimensional quantity). 
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2.1. Computational Method. To investigate the transition from free to forced 
convective flow, we developed a highly efficient solver for the governing equations. 
The numerical solution of forced convection problems becomes computationally 
very expensive due to the large linear systems which need to be solved for the 
velocity field. We have overcome this difficulty by using an eigenfunction expansion 
in the vertical direction using a Fourier sine series. 

The equations to be solved are Eqns.(5][7) subject to the initial and bound- 



ary conditions mentioned above. The velocity field is solved obtained using a 
streamfunction formulation, where the streamfunction ip is defined using the re- 
lation u = V X ipez- Applying the curl operator to Eqn. ([5|, we get the following 
Poisson equation for the velocity field: 

(15) 

The boundary conditions are given by "0 = at y = 0, ■0 = uq ^-t y = Ij V' = '^J-oU 
along the left boundary and dip/dx = along the outflow boundary. We now 
expand the streamfunction in a Fourier sine-expansion in the vertical direction as: 

N 

ip{x,y) ^ uoy + ^i^n{x)sin{mTy) (16) 



Substituting this expansion in Eqn. ( |15[ ) followed by a multiplication with sinirmry) 
and integration in the vertical direction, we get: 

^™ n^TT^tprn — ~2 ( sin(m'Ky)dy (17) 



dx^ Jq dx 

This solution method depends crucially on a theorem about Fourier transforms 
which guarantees convergence of the term-by-term integral of a Fourier series, even 



if the original series is divergent (see Theorem 9.8.1 from Prosperetti 2011 ). This 

dC 

dx 



issue is important to consider since the integrand ^ may not necessarily vanish at 



the boundaries. 

Eqn. ( |17[ ) is now a set of Ny ordinary differential equations in the horizontal 
direction, where Ny is the number of vertical grid-points. These equations can be 
discretised using the standard second-order finite difference method, resulting in a 
set of tridiagonal equations that can be solved efficiently. Once the velocity field is 
obtained from the streamfunction, the concentration field is updated using a second 
order explicit Adams-Bashforth method. During the very first time step, we use a 
second order Runge-Kutte method for advancing the concentration solution. 

The simulations are performed using 128 uniformly spaced nodes in the vertical 
direction and L x 128 uniformly spaced nodes in the horizontal direction. The vari- 
ables are located on the grid in a staggered fashion, with concentrations evaluated 
at the center of the cells and velocities evaluated on the faces. The streamfunction is 
defined at the top-right corner of each cell. We use constant intervals for timestep- 
ping, with a nondimensional timestep of At — 10"^, chosen to satisfy both the 
diffusive s tability condition < ^ and the Courant-Friedrichs-Levy condition 



< 1 (|Ferziger and Peric 2001 ). 



3. Results 



We have used the computational technique described above to investigate the 
transition from free to forced convection for Ra = 50, 100, 200 and 400, and for uq 
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Figure 2. Comparison between the theoretically predicted con- 
centration profile C{x,y) = 1 — erf{y^/Rauo/2y/x) (black line) 
and that computed from a simulation (blue symbols)) with Ra = 
400, uq — 10 along a line at the center of the CO2 source. 



varying between 0.1 and 10.0. The length of the CO2 source was varied between 
i = 10 to L = 40. In Fig. (|3|, we plot the Nusselt numbers computed for the 
simulations scaled with RaL, where the Nusselt number is defined as: 

Nu^ / - — iy = 0)dx (18) 

We first verify that the grid spacing used in the simulations is fine enough to 
resolve the boundary layer by performing a simulation with Ra = 400 and Mq = 10 ( 
the highest values for both parameters considered in this manuscript). We compare 
the profile of dissolved carbon dioxide in the vertical direction at the center of the 
plate with the similarity solution given by Eqn. (11). In Fig. ([2]), we plot the 



computed boundary layer profile using blue symbols and the analytical solution 
using the black line. The simulations are seen to match the theoretical solution well, 
confirming that the resolution is fine enough to resolve the boundary phenomena. 

From the scaling arguments presented in the previous section, Nu/RaL = 1/40 
in the free convection regime for low Rayleigh numbers and Nu/RaL = a/mq / RaL 
in the forced convective regime. These scalings are plotted using dashed lines in 
Fig. ([3]). It must be noted that the parameter combination RaL that occurs in 
the analysis is a Rayleigh-number that is based on choosing the size of the CO2 
source as a length scale instead of the depth of the medium. The simulations are 
performed upto time T = 10 and the Nusselt number is computed by averaging the 
instantaneous values from t = 9.5 to t = 10, and these average values are plotted 
using blue circles. It can be seen that both the computed Nusselt numbers and the 
critical velocity agree very well with those predicted using theoretical arguments. 
For concreteness, we present contours of the dissolved concentration of CO2 for 
simulations fixed with Ra = 400 and L = 10 in Fig. Q for background flow 
velocities of uq — 1, 2, 3, 4, 5 and 10. The predicted critical velocity for transition to 
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Figure 3. Comparison of scaled Nusselt number, computed 
from simulations (blue circles) and the theoretical predictions for 
free and forced convective regimes (dashed lines) for different back- 
ground flow rates, scaled as 

forced convection is Uc — 2.5. The contours show instantaneoous snapshots of the 
system a,t t — 10. It is seen that for velocities uq > Uc, the flow resembles forced 
convective behavior closely. Also presented in Fig. ^ are the computed Nu/RaL 
at different times for this problem with uq = 2 (blue line) and uq — i (black line). 
The simulation with = 2 is slightly below the critical velocity and demonstrates 
free-convective behavior which can be identified using the distinct minima in the 
dissolution rate followed by an enhancement due to the effect of fingering. The 
case with uq = 3 on the other hand achieves a steady state dissolution rate and 
does not display the minima. The dashed lines represent average Nusselt numbers 
that were used in Fig. Q. The fluctuations in the Nusselt number are due to the 
fact that the simulations are continuously forced (once every 100 At), to produce 
fingering phenomena when that is the preferred state of the system. In the absence 
of continous forcing, the simulations can reach an artificial steady state even when 
fingering is the preferred state after the initial perturbations have been flushed out 
of the domain by the background flow. 



The results presented in this manuscript describe how the presence of a back- 
ground flow alters the dynamics of CO2 dissolution in saline aquifers. The main 
implications of these results are: 

(i) An accurate understanding of the magnitude and direction of background 
flows is needed to locate monitoring wells. In the presence of a reasonably 
strong background flow, monitoring wells located right under the CO2 
plume to investigate dissolved carbon content might never record fluid 
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Figure 4. Contours of concentration of dissolved CO2 for simu- 
lations with Ra = 400 and L = 10. The predicted critical velocity 
is Uf. = 2.5, beyond which the dynamics near the CO2 source are 
predicted to resemble forced convection. The location of the CO2 
source is marked in red. 
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Figure 5. Nusselt number at different times for two simulations 
with Ra = 400 and uq = 2.0 (blue line) and uq — 3.0 (black line). 
The critical velocity from theory is uq = 2.5. 




with meaningful concentrations of carbon dioxide. In the presence of ex- 
perimental data on existing background flow conditions, these simulations 
can be used to suggest suitable locations to place monitoring wells, 
(ii) For low Rayleigh number systems, forced convection can result in far higher 
dissolution rates than free convection. Further, the magnitude of the rates 
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of dissolution can be estimated quite accurately using the well-known the- 
oretical formulae. Another benefit of the constant dissolution rate due to 
forced convection has to do with upscaling into larger reservoir simulators 
where the boundary layer cannot be resolved, 
(iii) In this work, we have not explicitly included the role of either hydrody- 
namic dispersion, or anisotropy of the permeability tensor. In the context 
of the simple theoretical arguments presented in this paper, hydrodynamic 
dispersion essentially serves to define an effective (higher) diffusion coef- 
ficient which will push the critical velocity lower and the Nusselt number 
under conditions of forced convection higher. Similarly, the presence of a 
permeability anisotropy with higher horizontal permeability will tend to 
make forced convection even more efficient than suggested in this manu- 
script. 
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